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Abstract 



<^ ' Vortex nucleation in a Bose-Einstein condensate subject to a stirring potential is studied nu- 

merically using the zero-temperature, two-dimensional Gross-Pitaevskii equation. In the case of 
: a rotating, slightly anisotropic harmonic potential, the numerical results reproduce experimental 

findings, thereby showing that finite temperatures are not necessary for vortex excitation below 
the quadrupole frequency. In the case of a condensate subject to stirring by a narrow rotating po- 



tential, the process of vortex excitation is described by a classical model that treats the multitude 
q \ of vortices created by the stirrer as a continuously distributed vorticity at the center of the cloud, 

but retains a potential flow pattern at large distances from the center. 
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I. INTRODUCTION 



Since the first experimental creation of a quantized vortex in a gaseous Bose-Einstein 
condensate (BEC) Q,Q], there has been significant experimental and theoretical advance in 
this field of research. However, even in these relatively simple systems, the understanding of 
such a basic issue as the mechanism for nucleation of vortices has not been straightforward. 

Several experimental methods of vortex creation are currently in use, including phase im- 
printing , cooling of a rotating normal gas 4], and conversion of spin angular momentum 
into orbital angular momentum by reversal of the magnetic bias field in a Ioffe-Pritchard 
trap QflQ. The topic of this paper is the rnechanical stirring of a Bose-Einstein condensed 
cloud with the use of optical and magnetic fields. This method, in turn, falls into two cat- 
egories: rotation of an anisotropic trap and stirring with a narrow potential, and we shall 
examine each of those in turn. 

References reported on vortex creation in gases that are confined in traps whose 

equipotential curves in the x-y plane have a slightly elliptic shape and are varied in time 
so that the axis of ellipticity rotates around the third axis. The formation of vortices 
in such traps is understood to be a consequence of the excitation of unstable quadrupole 
modes 0, Q] • Accordingly, vortices are seen to occur first when the angular frequency of 
rotation is in the vicinity of l/y/2 times the average radial trapping frequency lo (that is, 
;he frequency of the quadrupole mode divided by its angular-momentum quantum number) 



111 ]. However, the range of parameter values within which vortices are experimentally 



seen to be created is not yet entirely accounted for by theory. This shall be discussed in 
detail later. 

Finally, vortices can be excited by the use of a repulsive optical potential which is confined 
to a small region of the plane, thus mimicking a stick which is moved through the gas cloud 
in order to stir up vortices |l2j]. In this kind of geometry, local turbulence including the 
creation of vortex-antivortex pairs at the position of the stirrer is instrumental in inducing 
angular momentum in the gas [13] . 

In this paper, we set out to investigate to what extent the experimental findings can be 
accounted for by two-dimensional zero-temperature theory. We shall do so with means of 
numerical calculations and, when possible, classical hydrodynamical modelling. The paper 
is organized as follows. The equation of motion, the numerical methods and the elements of 
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the theoretical understanding of vortex nucleation are introduced in Section |H] In Sec. IIHI 
we study the formation of vortices in a rotating elliptical trap and compare our numerics 
with the experimental findings. Section HVI is concerned with the method of stirring the gas 
with a moving narrow optical potential. In Sec. we summarize and conclude. 



II. EQUATIONS AND NUMERICAL METHOD 

We assume the gas cloud to be at sufficiently low temperature and low density that it can 
be described as a pure condensate. Such a system is entirely governed by a scalar condensate 
wave function ^(r, t) that evolves in time according to the Gross-Pitaevskii equation 14. l^: 



The local fluid density and velocity are given in terms of the condensate wave function 
by n(r, t) = \ip(r, t)\ 2 and v(r, £) = (H/m) V(arg^(r, t)). The nonlinear term in Eq. ([T]) 
describes inter-particle interactions in the s-wave approximation and its strength is deter- 
mined by the s-wave scattering length a. V(r, t) represents the external potential and can 
in experiments be of both magnetic and optical origin. We shall in the following take V(r, t) 
to be the sum of two terms, a two-dimensional cylindrically symmetric harmonic-oscillator 
potential and an asymmetric time-dependent term: 

V(x, y, t) = l -mu)\x 2 + y 2 ) + AV(x, y, t). (2) 

The spatial extent of the one-particle ground state of the harmonic potential is called the 
oscillator length, a osc = (h/mto) 1 ^ 2 . Finally, the condensate wave function ip is in two 
dimensions normalized to the number of particles per unit length which we denote N/l z , so 
that l z is the average extent of the cloud in the z direction. 

We shall make frequent use of the Thomas-Fermi approximation, which neglects the 
kinetic term in the Gross-Pitaevskii equation ^(|. This approximation gives an accurate 
description of the bulk properties of the cloud when the coupling is strong; however, it fails 
to correctly describe short-range density variations such as that at the cloud boundary and 
close to a vortex core. The stationary-state density in the Thomas-Fermi approximation, 
assuming an isotropic harmonic-oscillator potential, is 

n(r) = no f 1 - ^jp^ (3) 



in the region where the right-hand side is positive and zero elsewhere. The central density is 
given by tiq = \imj (47rfr 2 a); the chemical potential \x is determined by normalization. In the 



two-dimensional case considered here, \x = 2hu^J Na/l z . The Thomas- Fermi cloud radius is 
.Rtf = 2 {Na/l z ) l / 4 a osc . Of importance is also the healing length £, which determines the size 
of a vortex core; in the Thomas- Fermi approximation it is given by £ = (a sc/-RTF) flosc- The 
condition of strong coupling is fulfilled when R TF ^> £ or, equivalently, when the coupling 
parameter Na/l z is large. 

In the numerical study, the Gross-Pitaevskii equation (P) is propagated in real time in the 
nonrotating (laboratory) reference frame, subject to an explicitly time-dependent potential, 
using the split-step Fourier method. We have chosen not to use any phenomenological 
model for finite-temperature dissipation such as an imaginary component in the time step 



(cf. liiLabd) 



As already mentioned above, we have chosen to study only two-dimensional systems, 
knowing that the physics of vortex formation for oblate systems will be well described by 
such a model. In fact, as we shall see, the model seems to be able to describe with limited 
accuracy the essential physics of vortex formation also for a prolate system. 

During the temporal evolution, we have, in addition to the complex wave function itself, 
monitored the total angular momentum, the mean extension of the cloud along suitably 
chosen axes (corotating with the disturbance) and the total fluid circulation within a region 
close to the condensate center. All these quantities are measurable experimentally, and we 
shall see that each of them exhibits a clear signal in its time evolution when vortices enter 
the interior of the cloud. 

For a vortex to be created, it must be energetically favorable compared with a nonrotating 
state [^l, but in addition, there must be a dynamical instability of the cloud that permits 
vortex creation. Generally, an excitation can only be created due to a moving disturbance 
in the fluid when the velocity exceeds a critical one, given by the Landau criterion j^ : 

v c = min — , (4) 
1 q w 

where q denotes the wave numbers of all possible quasiparticle modes. For a trapped BEC 

subject to a rotating perturbation close to the condensate edge, the Landau analysis yields 

a critical angular velocity in terms of the surface modes P. I22I: 

Q c = min y , (5) 



where / are the angular momentum quantum numbers of the surface modes. The result for 
typical trap parameters is a critical linear velocity v c slightly smaller than the speed of sound 
in the center of the trap c 23]. As we shall see, this simple analysis is not sufficient for the 
two kinds of geometry studied in the present paper, but needs to be modified in different 
ways. 



III. VORTEX CREATION IN A ROTATING ANISOTROPIC POTENTIAL 

Refs. j3, describe experiments performed on a Bose-Einstein condensed cloud contained 
in a harmonic-oscillator potential that has a small ellipticity e in the x-y plane which rotates 
with an angular frequency fl The time dependent potential can be written 

V(x, y, t) = l -muj 2 [(1 + e)X(t) 2 + (1 - e)F(t) 2 ] , (6) 

where the coordinates along the corotating axes vary in time as 

X(t) = x cos Qt + y sin Qt, Y(t) = y cos fit — x sin fit. (7) 

The twofold symmetric potential can only excite modes of quadrupolar symmetry. Therefore 
the Landau criterion, Eq. (jHJ), must be modified so that the summation is restricted to 
quadrupolar modes, resulting in the criterion 

a = f , (8) 

where u>2 is the excitation frequency of the lowest-lying quadrupole mode. In the Thomas- 
Fermi limit it is known to have the value uj 2 = \/2uj and therefore the Landau-criterion 
critical frequency for this trap geometry is Q c = uj/\/2 ~ 0.71 to. However, in order for 
vortices to be formed, a necessary condition is not only that the quadrupole mode is excited, 
but also that it is unstable 0, This has been shown to occur when the ellipticity is 



larger than the critical value 



The onset of instability above the quadrupole frequency has been investigated in detail 
in Ref. In Ref. 0, it was found experimentally that vortices were indeed created 

only for e > e when Q > u/y/2. On the other hand, for driving frequencies slower than the 



quadrupole frequency, vortices could still be created if the ellipticity exceeded a ^-dependent 
value that appears to be well approximated by the straight line e = 0.71 — Q/u. This can 
be thought of finite width of the resonance at Q = u/V2. In Ref. Q it was speculated 
that this finite width is either due to finite temperature or the excitation of surface modes 
with higher multipolarity m > 2. The former proposition seems unlikely because there was 
little temperature variation of the position of the lower threshold value. In addition, the 
experimentally found sub-resonance vortex nucleation is reproduced by the zero-temperature 
theory of Ref. j^J, as well as by the present numerics, as we shall see. Reference [24] puts up 
a pictorial model of vortex entry that focuses on the path taken by a vortex moving into a 
BEC cloud. The finite energy barrier that a vortex would have to overcome in order to move 
into a cylindrically symmetric cloud, even under external rotation, is found to be lowered if 
the cloud is elongated. The critical frequency Q is determined by the requirement that the 
energy at the saddle point which the vortex has to overcome is lower than the energy of a 
cylindrical, nonrotating cloud. The analysis thus applies to the case of a sudden switch-on 
of the rotation. The critical frequency is shifted down from the quadrupole frequency by 
an amount that depends on both e and the coupling strength. However, in the case of an 
adiabatically slow increase of the rotational force, the model of Ref. j^J would predict no 
vortex entry, because there will always be a local energy minimum for a vortex-free elongated 
state and therefore a finite energy barrier for a vortex to overcome. 

The numerical procedure mimics the experimental conditions of Ref. [sj], which has an 
effective two-dimensional coupling strength Na/l z = 56. The two-dimensional calculation 
neglects any degrees of freedom associated with the z direction, which is a valid approxi- 
mation since the experimental trap is oblate with an aspect ratio w/u z ~ 0.35 and in such 
a trap vortices are expected to be straight (25" 



|. The cloud is initially taken to be in 
equilibrium in a stationary trap with zero ellipticity in the x-y plane. The ellipticity e is 
ramped up from zero to its final value over a time t m i t = 80 uo^ 1 while the rotation frequency 
Q is held constant. The subsequent evolution is monitored for another 300 cu -1 ; typically 
the vortex nucleation starts after a time comparable to the ramping-up time when Q lies 
above the quadrupole frequency. To obtain a diagnostics for the presence of vortices, the 
ellipticity is subsequently ramped down between t = 300 c<j -1 and t = 380 uj^ 1 , whereafter 
the vorticity in an area enclosed by a circle of radius -Rtf is computed and averaged over a 
time span of 100cj _1 in the cylindrically symmetric trap. An average vorticity above 0.3 is 
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FIG. 1: Critical values of the anisotropy e of a rotating anisotropic trap for the onset of vortex 
nucleation. Asterisks denote the numerical results from solving the Gross-Pitaevskii equation, the 
dashed line is the analytical prediction of Ref . 9] , and the circles represent the experimental results 



ofRef. 



8]. The solid line represents the analytical estimate of Ref. [2], which is n I licory for sudden 



switch-on of the rotation, applied to the case of a chemical potential \i = Ylfiuj. 

taken to be indicative of the fact that vortices spend a nonnegligible time in the interior of 
the cloud; the results are not sensitive to the precise threshold value. 

Figure IT] shows the critical value of e as a function of Q. It is seen that the numerical 
reSU , tS are in g ood a g _ with t fie e„ t ai finding, and t fie t fie OTy „ t Ref. Q 
when Q is larger than the quadrupole frequency, but the agreement is somewhat poorer 
below. The qualitative features are, however, well captured by our numerics. This fact 
rules out the possibility that vortex excitation below the quadrupole frequency is purely a 
finite-temperature effect, since the present study is performed using zero-temperature theory. 



The prediction by Ref. 



24j is also in rough agreement with the experiment and the present 
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numerics but predicts a lower critical anisotropy, which is expected since it applies to the 
case of a sudden switch-on of the trap. We speculate that the present ramping-up time, 
although long compared to the period of the forcing Q -1 , is still not perfectly adiabatic. 
While the threshold values of e and Q agree with the experimental findings of Ref . jslj , the 
observed numbers of vortices turn out to be somewhat higher. This can be due to lack of 
equilibration or because our numerical counting method includes vortices on the edge of the 
condensate that are not visible to the eye. 

In order to better understand the transition to the vortex state, we study in detail the 
evolution of the z component of the angular momentum L z and the mean extension of the 



cloud along two axes corotating with the trap, X rms = y (X 2 ) and Y rms = y(Y 2 ). The 
deformation of the cloud is denoted S and defined as 

y 2 _ x 2 

a — ruus ' i ins (10) 



Y 2 + XL 



in accordance with Ref. 



24j ]. In Fig. |21 these four quantities, together with the vorticity 



N v in an area enclosed by a circle of radius Rtf are plotted against time for the choice of 
parameters Q = 0.78a; and e = 0.1. This point lies in the right part of the phase diagram of 
Fig. ^ above the quadrupole frequency; we shall soon comment on the difference in behavior 
for different Q. In this plot, the anisotropy e is initially ramped up for a time t = 80 a; -1 as 
before, but is then held constant throughout the simulation in order to display the physics 
more clearly. We see how the cloud is initially elongated and at the same time rotates with 
the trap. The cloud initially rotates 90 degrees out of phase with the potential, so that 
the elongation is in the direction of stronger trapping potential, consistent with Ref. 0. 
The elongation of the cloud would now begin to oscillate around an equilibrium value if 
the frequency and ellipticity were outside the critical region. Within the critical region, the 
maximum elongation instead grows larger until the outer edges of the cloud bend, close in on 
themselves, and capture vortices; the first vortex is seen to be detected just prior to t — 50a; -1 
when the cloud elongation is at its largest. After a turbulent period, the anisotropy of the 
cloud is somewhat decreased. During and after vortex capture, the orientation of the cloud 
is changed so that the elongation is in the direction of weaker trapping. The whole process 
is in accordance with previously reported results SB 0, Q . We do not observe an actual 
crystallization of the vortices into a lattice during the waiting time allowed for here. A recent 
preprint suggests that the crystallization of a lattice should occur even at zero temperature 
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FIG. 2: The upper panel shows the time evolution of the angular momentum (full curve) and 
number of vortices N v (dots) in a trap with an anisotropy e = 0.1 rotated at a frequency Q = 0.78 lo. 
The second panel from the top is a plot of the rms radii of the cloud along two axes corotating 
with the anisotropic trap: the full curve describes X ims = \/(X 2 ) and the dashed curve describes 
y rms = \J (Y 2 ). The third panel depicts the dimensionless cloud deformation 5 as a function of 
time. The four pictures on the bottom are density plots taken at the time instances t = 27a> -1 , 
t = 80a; -1 , t = 160a; -1 , and t = 400a; -1 , respectively. Brighter shades represent higher density. 
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FIG. 3: The upper panel shows the time evolution of the angular momentum (full curve) and 
number of vortices N v (dots) in a trap with an anisotropy e = 0.2 rotated at a frequency f2 = 0.64 to. 
The lower panel depicts the dimensionless cloud deformation 5 as a function of time. 

[rgj ]. but the difference in dimensionality, trap geometry and number of vortices prevents 
a quantitative comparison of the time scales. The effect of dissipative mechanisms on the 



lattice formation was studied in Refs. 



13 



3. 



Below the quadrupole frequency, the process is markedly different. Fig. |H] depicts the 
evolution of the angular momentum, vorticity, and cloud deformation 5 in this regime. The 
cloud is initially elongated, but the elongation of the cloud rotates in phase with the trap 
deformation so that 5 is positive. The cloud deformation is seen to grow to a maximum 
whereafter vortices enter the cloud, consistent with the physical picture of Ref. j^J. The 
buildup of angular momentum takes place gradually and smoothly compared with the situ- 
ation in Fig. 121 

Figure 0] displays the time evolution of the angular momentum for two different choices 
of Q, providing representative examples of the evolution of angular momentum below and 
above the critical line. Clearly, the physics of the vortex creation is different depending on 
whether the rotation frequency is larger or smaller than the quadrupole frequency. Above 
the quadrupole frequency (rightmost panel in Fig.EJ), the vortex formation is very fast, and 
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FIG. 4: Time evolution of the angular momentum in a rotating anisotropic trap for different 
values of the rotation frequency O, and trap anisotropy e. In the left panel, the rotation frequency 
ft = 0.64a; and the trap deformation e = 0.1 (full line) and e = 0.17 (dashed line). In the right 
panel, the rotation frequency O, = 0.78a; and the trap deformation e = 0.41 (full line) and e = 0.42 
(dashed line). For the solid lines, the parameter values are such that no vortices are created, while 
the oscillations in the case of the dashed lines are associated with vortex creation. 

the oscillations in the vortex-free regime are distinct from the turbulent behavior leading to 
vortex formation, although the values of e for the dashed and full curves differ by a very 
small amount. Below the quadrupole frequency, however, vortex formation is preceded by a 
gradual buildup of angular momentum, and the outgrowth of "arms" from the cloud seems 
to be absent. In this regime, it takes a finite time for the instability associated with vortex 
formation to occur. 

The critical lines can thus be determined with good accuracy above the quadrupole 
frequency, but below the quadrupole frequency, they have a more pronounced dependence 
on the waiting time. This fact together with the turbulent fluid flow and sensitivity to details 
that is connected with the vortex creation is the reason for the scatter of the numerical data 
points in Fig. ^ The gradual buildup of vorticity is consistent with the view that vortices 
enter the cloud only when there is an appreciable elongation (or equivalently, population of 
the quadrupole mode); this cannot be captured by the linear stability analysis of Ref. |lo| . 



An alternative, but probably equivalent picture, is the one proposed in Ref. 



2J, which still 



holds if the rotation is switched on moderately rapidly. Another possibility could be that the 
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vortex creation is dependent on the excitation of higher than quadrupolar modes due to trap 
imperfections (which in the simulations would then be due to the square numerical grid). But 
if this were the case, we would not have expected the fair quantitative agreement between 



(The occupation of modes during 
181 ]. but without special attention to 



numerics and experiment that is after all visible in Fig. 
vortex nucleation has been studied numerically in Ref. 
the difference between frequencies below and above the quadrupole frequency.) We therefore 
conclude that the cause for vortex excitation in this regime is not to be found in factors 
extrinsic to the zero-temperature Gross-Pitaevskii equation. 

IV. STIRRING A CONDENSATE WITH A LOCALIZED POTENTIAL 



In the experiment reported in Ref. |l2j, angular momentum was imparted to the con- 
densate by stirring it with a pair of narrow repulsive laser beams, placed opposite to each 
other and rotating about the center at a radius r s with an angular frequency Q, so that the 
linear velocity is v s = Qr s . The vortex formation is in this case governed by quite different 
physics compared to the previous section. In fact, there exists a crossover: in the limit of 
a pair of wide beams at large r s , the potential becomes again that of a rotating anisotropic 
trap and the analysis of the previous section applies. Here we shall, however, study the case 
of narrow beams at both large and small distances from the center. We take the potential 
to be a harmonic-oscillator one plus two symmetrically placed repulsive Gaussian "sticks" 
revolving about the origin: 

1 1 1, ( l r - r 0( t )| 2 |r+r (t)| 2 \ 

V(x, y) = -mu 2 {x 2 + y 2 ) + V U + e J . (11) 

The stirrer coordinates are given by r (t) = (xo(t),yo(t)) with 

xo(t) = r s cos fit, 

y Q (t) = r s sinfit, (12) 

and the stirrer width o is smaller than the cloud size Rtf but not smaller than the healing 
length £. We have chosen to consistently work with a pair of stirrers in order to facilitate 
comparison with the experimental data of Ref. Q]. We have, however, found no qualitative 
differences when comparing with test runs made with one stirrer. 

The physics of the vortex creation is the following Q|. When an object or a repulsive 
potential is moved through the condensate, pairs of vortices with unit positive and negative 
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circulation are created in the wake. Positively oriented vortices will be created to the left of 
the stirrer relative to its direction of movement, and negatively oriented vortices, antivortices, 
to the right. (This can be realized by considering how the fluid flows back around the stirrer 
to fill the wake.) The vortices and antivortices follow the local current which initially is 
directed at an oblique angle, away from the stirrer and in the direction of movement. The 
antivortices will quickly drift out of the system while the positive vortices may enter into the 
cloud depending on whether they have enough energy to do so (assuming that the stirrer 
rotates counterclockwise around the center; otherwise the role of vortices and antivortices is 
interchanged). The whole process is turbulent and is complicated by inter- vortex interaction, 
sound waves and density variations. 

The conditions for vortex creation are thus twofold: the linear velocity of the stirrer 
must be large enough to create vortices, and the angular velocity of the stirrer must be 
large enough to contain the vortices in the interior of the cloud. We study first the former 
requirement. Vortex pairs are only created when the stirrer velocity exceeds a critical value 
given by the Landau criterion, Eq. ©. For the situation in Ref. |l2|, the prediction of Ref. 
[23j is v c = 0.24 c. However, this analysis was carried out for a stirrer that moves close to 
the edge of the condensate. For the case of a Gaussian stirrer that moves in the bulk of the 
condensate, the abrupt potential gradient is expected to result in a lower critical velocity v c 
than that given by the surface- mode calculation [27]. In this case, a critical velocity about 
an order of magnitude smaller than the sound velocity c has been observed numerically [28 1 
and experimentally jl^ . I29I ] . An approach that inserts the calculated energy and momentum 
of a vortex pair into the Landau criterion, Eq. (j5J, has yielded a critical velocity of the 
correct magnitude 3(3]. 

We have numerically investigated the critical velocity for the case when the coupling 
parameter Na/a osc = 56, giving a Thomas- Fermi radius -Rtf = 5.5 a osc , and a stirrer moving 
at a radius r s = 0.2 .Rtf from the center. When the stirrer size a = 0.05 Rtf, the critical 
velocity v c ~ 0.45 c and when a = 0.1 -Rtf, v c ~ 0.35 c. Due to the noisy character of the 
data, the error bars on the observed critical velocities are about ±0.1c. This critical velocity 
is larger than that of Ref. [12], but the difference in coupling strength is quite drastic: due 
to numerical constraints we have chosen to operate with a coupling strength Na/l z that is 



smaller than that of Ref. 



12j by approximately a factor 10. 



When the stirrer velocity is well above the critical one so that the number of vortices 
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is large, the build-up of angular momentum resulting from the complex process of vortex 
creation and entry can be described classically. The torque exerted on the fluid by the 
moving disturbance is estimated as follows. Assuming that the stirring potential Vq is strong 
enough, the fluid that is instantaneously right in front of it is accelerated to the velocity of 
the stirrer v s . The momentum increase per unit time due to one stirrer is therefore equal to 
the product of the mass of the fluid occupying an area swept out by the stirrer in unit time 
and the increase in velocity of the fluid at radius r s : 

dp 

— = 2av s mn 2D {r s ){v s - v(r s )), (13) 

where U2d{t) is the two-dimensional density of the fluid; ti2d(i") — n&Dirflz- The width of 
the stirrer is taken to be 2a. The momentum increase is azimuthally directed, and therefore 
the torque acting on the gas is r s times the quantity on the right-hand side of Eq. (fT^|) . 
When two stirrers are used, the momentum increase is doubled and the total torque is 

t = 2- 2av s r s m(v s - u(r s ))n 2 o(r s ). (14) 

The effect of the ensuing turbulent behavior, including the creation of vortex pairs, is to 
quickly distribute the momentum over the whole cloud, and when the number of vortices 
is large, the fluid rotates on average as a solid body: v (r) = lot. We can assume that the 
redistribution of momentum happens instantaneously, because the characteristic timescale 
for the spreading of the disturbances is t spre ad = -Rtf / c. For realistic parameter values, this 
timescale is indeed much shorter than the timescale for the growth of the angular momentum, 
and therefore the system can at all times be assumed to be in solid-body rotation. However, 
during the stirring, we do not expect the vortex array to extend outside the radius r s of the 
stirrer. Outside the region enclosed by the path of the stirrer, the fluid velocity is the sum of 
velocity fields from the individual quantized vortices, decaying as the inverse distance from 
the vortex. A reasonable assumption for the instantaneous fluid velocity is therefore 

\V (15) 

where the instantaneous solid-body rotation frequency u(t) characterizes the rotation rate. 
Multiplying the local velocity, Eq. ()15j) . with r and the Thomas- Fermi density, Eq. ([3*]). and 
integrating over the whole cloud, one obtains the total angular momentum per particle at 
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time t, 

L(t) = mu{t)r 2 s F j , (16) 

where the function F(x) = 1 — x 2 + x 4 /3. The derivative of the angular momentum is to be 
equated to the torque, Eq. (|14jl. and the substitution of the angular momentum L for the 
fluid velocity v(r s ) according to Eq. (fTT)|) results in 

L = t?( m \ {™r s v s F{r s /R TF ) - L) . (17) 
b {r s /KTF> 

Solving Eq. (fTTj) for L(t), one obtains the time dependence 

L{t) = L fin (l - e- 2Rt ) , (18) 
where the final angular velocity is 

^ = ™Hi£) > + 5 (£)')■ (19) 

and the rate constant: 

2av s n 2 D(r s ) 

= F(r s /R TF ) ■ (20) 
The factor 2 multiplying the rate constant in Eq. ()18|) is again the factor that reflects the 
number of stirrers. The torque is given by 

r = 2 ■ 2mar s v 2 s n (2D) {r s )e~ 2m . (21) 

The maximum value of the final angular momentum Lg n always occurs when r s /i?TF = 



(9 — v21)/10 ~ 0.66. The fact that the angular momentum is optimized when the stirrer 
is halfway between the center and the edge of the cloud is simply due to the fact that the 
linear velocity of the stirrer has been held constant. The maximum angular velocity is thus 
determined by a balance between on the one hand the increase in the amount of fluid that is 
enclosed by the path of the stirrer when r s is increased, and on the other hand the decrease 
in the angular velocity Q. If Q is held constant instead of v s , then the angular momentum 
will indeed be maximized by putting the stirrer at the edge of the condensate, as can be 
seen by substituting Qr s for v s in Eq. (fTT)|) and differentiating. 

The classical model is valid in the case when many vortices are created, so that the 
assumption of solid-body rotation holds. The predicted final angular momentum L fin , Eq. 
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(|19|). should therefore be much larger than unity; since the stirrer radius r s is on the order 
of the Thomas- Fermi radius -Rtf, this translates to the criterion 

' I ^ 1/4 



» tH • ( 22 ) 
a osc uj \NaJ 

The classical model is thus easily satisfied when the coupling is strong, as expected. In 
addition, the model is expected to break down when the stirrer is in the vicinity of the 
cloud radius, r s m -Rtf, because of the Thomas- Fermi approximation. The rate R is also ill 
described by the classical model when it is small; requiring R 3> t~^ read (see below Eq. (JHJ)), 
we obtain the condition 

>rV77T- (23) 



-O'oscJ Va osc u;/ l z V Na 
The results of the numerical calculations are found to agree well with the classical model. 

Calculations have been performed on a two-dimensional system with Na/l z = 40, yielding 

a ratio between the Thomas- Fermi radius and the healing length -Rtf / £ ~ 25, so that 

the Thomas- Fermi approximation is applicable. The sound velocity is c ~ 3.6co>a sc- The 

stirrer velocity has been varied between u>a osc and 3ooa osc . Fig. displays the final angular 

momentum Lg n of a two-dimensional cloud as a function of r s with fixed stirrer velocity 

v s , as well as the initial growth rate .R (found by fitting the data to a function of the form 

of Eq. (|18jl). An illustrative example of the variation of angular momentum with time for 

a particular value of r s is displayed in Figure |H1 We have chosen the parameters so that 

the agreement is fairly good, but not perfect. However, the qualitative appearance of the 

L curve is the same for all parameter values. The data is noisy because of the turbulent 

process, especially for the rate constant which is difficult to determine accurately, but on 

average it is seen to agree fairly well with Eqs. (j!8H20j) for higher values of v s and a so that 

the conditions (|22H23j) are met. 



The classical model is seen to describe well the experiment reported in Ref. [12j], even 
though that experiment was performed in a prolate trap with anisotropy u/u z — 4.3. For 
this trap, the Thomas-Fermi radius in the radial direction is .Rtf = 12.6a osc , and the speed 
of sound in the center of the trap is c = 8.5u;a osc , where u and a osc are the radial trap 
frequency and oscillator length, respectively. With a linear stirrer velocity v s = 0.2 c, the 
number of vortices N v was recorded as a function of the stirrer radius r s . In the limit of a 
large vortex density, the fluid is on the average in solid-body rotation, for which L = HN v /3. 
However, we are here operating with a relatively small number of vortices, the maximum 
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being about 20 vortices in the experiment. Corrections to the solid-body result may therefore 
be important. Indeed, by calculating the angular momentum for a triangular vortex lattice 
in a cylindrically symmetric Thomas- Fermi cloud for different values of the lattice constant, 
we find the corrected number of vortices as a function of angular momentum to be well fit 
by the linear relation 

N v « 3- - 6. (24) 

This formula obviously breaks down when the number of vortices is close to 1, but is accurate 
for N v > 10. Inserting the numerical values into Eq. ([T9j) and using Eq. (|2*3Jl for the number 
of vortices, we obtain for the specific conditions of Ref. 3| the relation 

N v = 6Ax (l-x 2 + ^x 4 ) - 6, (25) 

where x = r s /R^-p- The peak value of this function (at x = 0.66) is approximately 20.5. 
This result is displayed in Fig. and it is seen that the agreement is good, despite the 
drastic difference in geometry (prolate and two-dimensional, respectively). 



V. CONCLUSIONS 



The nucleation of vortices in Bose-Einstein condensed gases contained in time-dependent 
traps has been studied numerically for two different kinds of geometry. In the case of a 
rotating, deformed harmonic potential, vortices are nucleated when the anisotropy of the 
trap exceeds a critical value. Above the critical frequency for excitation of quadrupole modes, 
the numerically calculated critical anisotropy is in accordance with previous theoretical and 
experimental results. Below the quadrupole frequency, the numerical results are again in 
accordance with previous experimental results, showing that the nucleation of vortices in 
this regime can indeed take place at zero temperature. For a cloud stirred by a narrow 
moving repulsive potential, a classical model has been found to describe accurately the 
process of angular-momentum buildup. The classical argument is based on the fact that the 
velocity pattern is approximately a solid-body rotation in the interior of the system, where 
the density of vortices is high, but potential outside the path of the stirrer. As a result, the 
angular momentum approaches its final value exponentially. This model is applicable when 
the number of vortices is large, i. e. when the stirrer velocity and width are not too small. 
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FIG. 5: Final angular momentum (upper panels) and initial rate of increase of angular momentum 
(lower panels) for a condensate stirred by two rotating laser beams, as functions of the distance of 
the stirrers from the trap center. The stirrer velocity v s , stirrer width a and peak stirrer potential 
Vo are stated above each panel in units of the trap frequency oj, oscillator length a osc , and trap 
energy too. The lines are the classical equations explained in the text and the asterisks connected 
with lines are numerical results from the two-dimensional Gross-Pitaevskii equation. The coupling 
strength is chosen to Na/l z = 40, whereby the sound velocity c ~ 3.6 and Thomas-Fermi radius 
Rtf ~ 5.0 in trap units. The agreement with the classical prediction is best in the panels on the 
lower right, where the number of vortices is the largest. 
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FIG. 6: Time evolution of the angular momentum and the number of vortices of a condensate 
stirred by two rotating laser beams. The solid line denotes the numerically computed angular 
momentum, the dashed line is the classical prediction, the dotted line is a fit of the solid curve 
to a function of the form (|18|) . and the dots signify the number of vortices in the interior of the 
cloud. The parameter values are the same as in the lower left two panels of Fig. \E\ v s = 3uia osc , 
a = 0.4a osc , and we have chosen r s = 0.6-Rtf- 
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FIG. 7: Number of vortices after stirring as a function of the stirrer position for the conditions of 
the MIT experiment of Ref. The solid line depicts the classical model and the asterisks with 
error bars are the experimental results. 
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